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■ We analyzed the stochastic behavior of systems controlled by autocatalytic reaction A + X — » 

X + X. Assuming the distribution of reacting particles in the system volume to be uniform, we 
introduced the notion of the point model of reaction kinetics, and derived a system of differential 
£SJ . equations for probabilities of finding n = 0, 1, . . . autocatalytic particles at a given time moment. It 

, | has been found that the kinetic law of the mass action cannot be supported by stochastic model. 



< 



CO 



O 
o 



> 

O 



e 

T3 



PACS numbers: 02.50.Ey, 05.45.-a, 82.20.-w 
Keywords: stochastic processes; autocatalytic reactions; lifetime 



' Introduction 

o 

D 

The system denoted by S is defined as an aggregation of particles capable for autocatalytic reactions. Symbols X 
and A are used for notations of autocatalytic and substrate particles, respectively. The system is "open" for particles 
A, when it is in contact with a large reservoir which keeps constant the number of particles A. In contrary, the system 
is "closed", when particles cannot be injected in or extracted from the system. We suppose that the system is always 
closed for the autocatalytic particles X. 

As known, there are many papers and books 0, 0, EJ Q dealing with the influence of spacial (mostly diffusive) 
motion of particles on the kinetics of reactions, however, in many cases the complexity of the problem did not allow to 
obtain exact results. In order to analyze one of the decisive factors determining the stochastic nature of autocatalytic 
""O 1 reactions, in the present paper we do not deal with the diffusion of particles, instead we introduce the notion of the 
point model of reaction kinetics, similarly to that used in the theory of neutron chain reactors. It is to note that 
the point model of the reaction kinetics is based on a sever assumption according to that the reacting particles are 
distributed uniformly in the whole volume of the system. In this case, it can be stated that the probability of a 
reaction between two particles is proportional to the product of their actual numbers in the system. This approach 
is a highly simplified but useful model for systems, in which the number of particles is small, and the largest distant 
within the system is smaller than the characteristic length of interaction between the particles. 

I the sequel we will study closed systems controlled by random process in which the particles X catalyze the 
conversion of particles A into further X 's. Changing adiabatically the number of particles A with help of a suitable 
■^j- , reservoir, we can investigate the influence of the environment on the system behavior. 

We say that a system is in a living state when the conditions for autocatalytic reactions are existing, while a system 
is in a dead state when there is no possibility for autocatalytic reactions. The random time spent by a system in 
the living state is the lifetime of the system. One of the aims of this paper is the determination of the probability 
distribution of the lifetime and the study of its random properties. 

The organization of the paper as follows. To elucidate the approach that we will follow in the present work, in 
Section 1 we deal with systems controlled by the autocatalytic reaction A + X — > X + X. We derive and solve a 
system of differential equations determining the probabilities to find n = 0,1,... , Na new X particles arisen from 
A particles during the time interval (0,i). Since the system is closed, the number of particles A is decreasing with 
time, and finally becomes zero, i.e. the system reaches its dead state. By using the generating function equations, 
O ■ in Section 2 we show that the stochastic approach brings about an equation for the mean value of the number of 
autocatalytic particles completely different from the classical rate equation based on the kinetic law of the mass action. 
*h , In order to obtain some kind of solution of the hierarchical equation system derived for moments, the moment- closure 
(decoupling) approximation has been applied by many authors however, the consequences of the this procedure 

were hardly investigated. Since we succeeded to obtain exact results for the time dependence of moments, in Section 3 
we analyzed the error caused by the decoupling approximation. Finally, in Section 4 we study the random properties 
of the system lifetime. 
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I. SYSTEMS CONTROLLED BY REACTION A + X -> 2X 

We would like to study the properties of closed systems S of constant volume V controlled by the reaction 

A + X h X + X. 

Denote by N A and Nx the numbers of particles A and X, respectively at time instant t = 0. As the particles X 
convert A into further A's, so the number of particles A strictly decreases, and finally, the autocatalytic process is 
stopped. Obviously, in this case the system contains only X particles, and this is the state which is called dead state. 

A. Kinetic rate equation 

Denote by m(t) the expectation value of the number of particles A converted into X during the time interval [0, t]. 
Obviously, 

m A (t) = N A - m(t), and m x (t) = N x + m{t) (1) 
are the average number of particles A and X, respectively at the time moment t > 0. Introducing the notations 

M m (*) fn\ Na a (n\ Nx 

a(t) = -^, c A (0) = — , and Cx {0) = —, 

according to the kinetic law of mass action we can write the rate equation in the following form: 

^- = k c [c A (0)-a(t)} [c x (0) + a(t)}. 

Taking into account the initial condition a(0) = 0, and using the abbreviation u = k c t/V, one obtains that 

1 - e^ " 
1 + pe 

where 



a(u) = c A (0) - ; ___ Nou , (2) 



and 



N A 

P=— and N = N A + N X . 
i\ x 

For the sake of comparison with the results of calculations of the stochastic model, let us introduce the ratios: 

m A (u) e^"" 
rA(u) = — = (l + p) l + pe - NaU > (3) 

m x (u) 1 
rx{u) = — = {l + p) l + P e-^ - (4) 

B. Stochastic model 

Denote the number of particles A converted into particles X during the time interval [0, t] by the random function 
e Z+, where Z + is the set of nonnegative integers. At the same time, note that is the number of new X 
particles appearing in [0,t\. Let 

V{li(t)=n\Z(t)=0}=P(t,n), 0<n<N A (5) 

be the probability to find n new X particles at the time moment t > in the system 5 provided that at t = the 
number of new particles X was zero and the system contained N A particles of type A and N x particles of type X. 
The pair of integers {N A , N x } denotes the initial state of the system S. Assume that a At + o(At) is the probability 
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that any of particles X brings about a reaction with any of particles A. If the numbers of A and X particles in the 
system at time instant t > are N A — n and Nx + n, respectively, then 

h n At + o(At) =a[N A - n] [N x + n]At + o(At) (6) 

is the probability that the reaction A + X 2 X is realized in the time interval (t, t + At). It is to note that 
a = k c /V and Hna = 0- By using the expression © one can write 

p(t + At, n) = p(t,n)(l - h n At) + h n -ip(t,n - l)At + o(At), 

and obtain immediately the equation 

dp(t, n) 
dt 

If n — 0. then 



-h n p(t, n) + hn-ip(t,ri— 1), 1 < n < N A . (7) 



— — — = -h o p{t,0), (8) 

and taking into account the initial condition p(0, 0) = 1 one has 

p(t,0) = e- hot = c - NaNx u . (9) 

By introducing the Laplace transforms 



U n {s)= / e- st p(t,n) dt, n = l,2,...,N A , 
Jo 

since p(0,n) = 0, if n > 0, we obtain from Eq. Q the recursive relation 

(8 + h n )U n = hn-iU n -i(s), n>0, 
which can be easily solved by starting with 



0b 00 = -^t-- (io) 

s + h 

The solution can be written into the form: 

Un(s)= , — ^r^7 , hy n=l,2,...,N A . (11) 

In order to obtain the probabilities p(t, n), 1 < n < -/V^, the next step seems to be very simple: one has to expand 
U n (s) into the series of partial fractions, and then to perform the inverse Laplace transformation. However, the task 
is slightly more complex. 

If N A — Nx = n c < 0, then it is obvious that the quantities ho, hi, . . . , h n are all different and so, we can write 



n „ 

^HEixb (12) 



where 



(h - h k )(hi - hk)- ■ ■ (hk-i - hk){hk+i - h k ) ■ ■ ■ (h n - h k ) ' 
n = 1, 2, . . . , N A and k = 0, 1, . . . , n, 

and since Coo = 1, we obtain 

n 

p{t,n)=Y,C nk e- h " t , n = 0,l,...,N A . (14) 

k=0 
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If Na — Nx — n c > 0, then it can be easily proved the following statement: if n c is even, i.e. n c = 2k c , then 

h nc — ho, h nc -i = hi, . . . , hk c +i = hk c -i, 
while if n c is odd, i.e. n c = 2k c + 1, then 

h nc = 7 h Uc -i = hi, ... , /ife c +i = hk c ■ 
Let us denote the denominator of U n {s) by 

V n (s) = (s + h )(s + hi)---(s + h n ), 

and introduce the notation 
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Q{j,s) = \{( s + h l f. 



V n {s) 



i=0 

In the case of n > n c we can show that 

Q(k c - 1, s)(s + h kc )(s + h nc+ i) ■ ■ ■ (s + h n ), if n c = 2fc c , 
Q(k c ,s)(s + h nc+1 ) ■ ■ ■ (s + h n ), if n c = 2k c + 1, 

while in the case of n — n c we obtain 

Q(k c - 1, s)(s + h kc ), if n c = 2k c , 
Q(k c ,s), ifn c = 2fc c + l. 

When n < n c and < n < k Cl then we have 



V ne (s) 




FIG. 1: Probabilities of finding n = 0, 1, 2, 4, 5 new X particles at time moment t provided that the initial state of the system was {5,1}. 
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Vn(s) = Y[{S + hi) 
i=0 

for both, even and odd values of n c , but in the case of k c < n < n c we should use the formula: 

n™=o™ _1 ( s + hi) Q(n -kc — l, s)(s + hk c ), if n c = 2k c , 
nr=o" _1 ( s + h i) Q( n - fcc - 1. s), if n c = 2k c + 1. 



V n (s) 



For the sake of illustration it seems to be useful to write down the equations of U n (s) when = 5 and Nx = 1- 
One obtains 

TT , \ hlh\h 2 , , 1 



(s + h ) 2 (s + h 1 ) 2 (s + 112)3' uw s + h 

TT I \ h^h^hl TT { \ ^0 



{s + h ) 2 {s + h 1 ) 2 (s + h 2 y w {s + h )(s + hx)' 

TT i \ hoh\h 2 . . hohi 



(s + h )(s + h 1 ) 2 (s + h 2 )' zw (« + ho)(« + /ii)(« + h3)" 

One can see in FlC^that the probability p(t, n) versus t curves with exception of n = and n = 5 have a well-defined 
maximum the location of which is increasing with t. 

II. GENERATING FUNCTION EQUATION 

In order to make easier the derivation of equations for moments of the random function £(t) it is worthwhile to 
introduce the generating function 

N A 

g(s,z) = J2Un(s)z n . (15) 

n=Q 



It can be proved that g(s, z) satisfies the following equation: 

az 2 (l - z) d - a(n c - l)z(l - z) ^£lfl - [ s + h (l - *)] £(s, z) + 1 = 0. (16) 



In many cases the use of the exponential generating function 

N A 



g exp (s,y) = Y,Un(s)e n y (17) 



(18) 



n=0 

is more advantageous, therefore we write it down also: 

_ d 2 g ■ (s,y) _ dg ( s ,y) y) - 1 = 0. 

ay ay 

We can see that 

lim g(s, z) = lim g exp {s, y) = -, 
*->l y^O " s 

and it is elementary to show that the generating function 

g(t,z) = £^ 1 {g(s,z)}, 

where C^ 1 is the operator of the inverse Laplace transformation, satisfies the equation 

9g(t,z) ,dg(t,z) 2 d 2 g{t,z) 

— -Q t — = -fto(l - z) g{t, z) - a(n c - l)z{l - z) — ^- V az (1 - z) — — . (19) 
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with 

g(0,z) = l, and g(t, 1) = 1. 
Similarly, the exponential generating function 

9ex P {t, y) = £~ 1 {g e x P (s, z)} 
is nothing else then the solution of the equation 

= _, o(1 _ eV) gexp{t , y) _ anc{1 _ e .) ^SffiM + a( i - e «) (20) 

with 

gex P {0, y) = l, and g eX p(t,Q) = 1. 



A. Expectation value 

The numbers of particles X and ^4 at time instant f > are given by 

U(t) = N x +£(t), and U(t)=N A ~at), 
respectively. From these one obtains 



where 



E{£(t)}=mi(t) = 



mx(i) = 


-- N x + 




mx(t) = 


-- N A - 


mm, 


~dg(t,z)~ 




'dg e xp(t,y)~ 


dz 


2 = 1 


dy 



(21) 
(22) 



(23) 



By using the equation (|20|l after elementary algebra we obtain 

dm\(t) 



dt 

where 



ho + an c mi(t) — am^t), (24) 



m 2 (t)^E{C 2 (t)}, 

and this equation can be rewritten in the form: 

dmi(t) 



dt 



= a[N A - mi(t)] [JV X + rrn(t)] - aD 2 {£(t)}. (25) 



As seen the appearance of the variance D 2 {£(i)} brings about the loss of validity of the kinetic law of the mass action. 
If D 2 {£(£)} w 0, i.e. if m^if] ~ [ttzi (i)] 2 , then the deterministic kinetic equation can be applied. 

The equation (|24() shows the hierarchical structure of momentum equations. If we would like to calculate m\(t) by 
solving the equation (|24|) . we see that for this we need the second moment m,2{t). The equation of the second moment 
contains the third moment etc. In this way we can obtain the time dependence of the average number of X particles 
mi(t) when we solve the hierarchical system of equations which seems to be rather difficult. Therefore, in the practice 
the method of decoupling 1 has been often used for finding an " approximate" solution, however the consequences of 
this procedure were hardly investigated. In the next Section we try to give some elementary analysis of the problem. 



1 The decoupling is the substitution of the moment mfc(t) with an expression containing moments with indices smaller than k. 
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FIG. 2: Expectation values of the relative numbers of particles of type X (upper curves) and of type A (lower curves) versus time provided 
that the initial state of the system was {5, 1}. 



In the present case it is fortunate that we can determine the probabilities p(t, n), n = 0, 1, . . . , Na from and 
so we can calculate directly the exact time dependence of moments 

N A 

m k (t) = E{[£(t)] fc } = n k p(t, n). (26) 

n=0 

The expectation values of the relative numbers of X and A particles, i.e. the mean values of the random functions 

px(t) = l + fg, and PA(t) = l-^± 

versus time can be seen in FIG. [21 The upper curves refer to the X, while the lower curves to the A particles. The 
continuous and dashed curves show the results of calculations for the stochastic and deterministic models, respectively. 
As expected the difference between the stochastic and deterministic description is relatively small at the beginning 
and at the end of the process. 



B. Variance 



We have seen that the occurrence of the variance D 2 {£(i)} in the differential equation of l|25|l is the source of the 
invalidity of the kinetic law of the mass action. Therefore, it seems to be worthwhile to look at the time dependence 
of the variance and relative variance of the number of X particles. 

Since £*•(*) = N x + and = N A - £(*), it is trivial that 

D 2 Ux«} = D 2 {a(i)} = D 2 U(i)}, 

where 

N A 

D 2 {^(t)} = ^[n-m 1 (t)] 2 p(t,n). 

n=0 



By using the inverse Laplace transform of the expression Qllfl. we can calculate the variance D 2 {£(£)}. The results 
of calculations are plotted in FIG. El Observe, the increase of the number of starting X particles from 1 to 3 
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FIG. 3: Time dependence of the variance D 2 {£x(t)} and the relative variance D 2 {£x (t)}/E 2 {£x (t)} in two initial states: So = 
{5,1}, {5,3}. 

brings about a significant decrease of the relative variance. The fluctuation of £(t) increases sharply at the beginning 
of the process, and decreases slowly with increasing time. After elapsing a sufficiently large time we see that the 
fluctuation is negligible, and it means that the kinetic rate equation can be regarded in long time limit as an acceptable 
approximation. 

III. MOMENT-CLOSURE APPROXIMATION 

In order to show the consequences of heuristic methods applied often for obtaining an " approximate" solution of 
the hierarchical equation systems, let us derive the exact equations for m\{t) and m-zif). It follows from (|2U|I that 



dmi {€) 
dt 



= ho + an c mi(t) — a m,2(t) 



and 



dm 2 (t) 
dt 



= ho + (2/io + an c ) mi(t) — a{2n c — 1) m 2 (t) — 2a m 3 (t). 



(27) 



(28) 



We see that further equations are needed for m^{t) 1 m^{t) 1 . . . , and this makes the problem rather complex, however, 
the stationary values = lirut-joo rrik(t), k — 1,2,... can be easily obtained. Since m[ st ^ = Na, from l|27|l it 

follows that m 2 st ^ — N A , and by induction it can be proved that 



lv Ai 



k = l,2,. 
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If we would like to determine the time dependence of the mean value mx(t), then we should solve the hierarchial 
system of moment equations. The usual treatment of this problem is to use the moment-closure approximation, 
i.e. to truncate the hierarchical equation system. However, this procedure may bring about hardly foreseeable 
consequences. 

For example, if we substitute 7713(7;) in Eq. I|28|) with the product m\(t) 7712(7;), then we obtain two differential 
equations for m\(t) and 7712(7;), namely 2 

dmi(t) 

— ho + an c roi(t) — a 777,2(7;), 

at 

and 

— ^ — ho + (2/io + an c ) 777,1(7;) + a(277 c — 1) 7772(7;) — 2a mi(t) 7772(7;). 
Eliminating 7772 (t) in the second equation we can write that 

du 2 + 2n c[N A - m x (u)] [N x + mi(u)] + (1 - 3n c ) ^ ' + 2roi(u) ; = 0, (29) 

where u = at. We have to solve this equation, by taking into account the initial conditions 

dmi(u) 



mi (0) = 0, and 



= N A N X . (30) 

u=0 



du 

For the sake of simplicity we assume that Na = Nx — M, i.e. n c = 0. In this case we obtain 

d ( dmi 2\ 

It follows from this 

dm 1 ^ 

mi + ni 1 = C, 



du 

where C can be determined by applying the conditions (|30|) . As seen, C = Af 2 , and so we have 

— - = - im\ + 777! - M 2 , 

du 

the solution of which is given by 

1 _ p— Am« 

mi(u) = 2M 2 - — ; 5 — ; — , (31) 

v ; Am (1 + e- x ^ u ) + 1 - e - x *' u y ' 



where Xm — V4M 2 + 1. The first thing that can be seen immediately is 

lim 7711(77) = M 



y/l + l/(2Af) 2 + 1/2M 

and this limit value does not correspond to the exact relation limn^oo 7771(77) = M. 

In the upper part of FIG.^lwe plotted the exact "ex" and approximate "cl" values of E{px(u)}) versus u — at, and 
for the sake of comparison the curve "rt" obtained from Eq.Q is also plotted. It is remarkable that the decoupling 
7773(7;) w mi(t)m2(t) brings about larger error than the rate equation Q which corresponds to the decoupling 7772 (t) « 
mi(t)mi(t). As seen, the difference between the curves "ex" and "cl" vs. u is negligible at the beginning and at the 
end of the process. The lower part of FIG. 01 shows the time dependence of the differences between the curves "ex" 



2 Though these equations are different from the exact equations 1271 and 1281 . for the sake of simplicity, we are using the same notations 
for the approximate first and second moments. 
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FIG. 4: In the upper part the time dependence of E{px(w)} calculated by exact "ex" and approximate "cl" methods is shown. For the 
sake of comparison the curve " rt" obtained from Eq. J3J is also plotted . In the lower part the difference between the curves " ex" and " cl" 
vs. u = at is shown for three values of M . 



and "cl" for three values of M. Note, the difference decreases sharply with increasing M. Summarizing, we can state 
that the decoupling m 3 (t) « mi(t) m^ijt) brings about an unacceptable error in the calculation of mi(t). 

However, there is an other way to truncate the hierarchical equation system. We can calculate the cumulants of 
£(t) from the derivatives of K(t,y) = log g exp (t,y) at y = 0. As known, the cumulants 



K k (t) = 



d k K(t,y) 
dy k 



are zero in the case of normal distribution if k > 3. Although the random function does not follow the normal 
distribution at any t > 0, in contrary to this we assume that n^{t) — 0, and hence we use the following decoupling: 



m 3 (t) w 3m 2 (t)m 1 (t) - 2m\(t). 



(32) 



This closure is called by some authors |6j "normal approximation". In order to simplify the calculations, we are 
dealing with the case when Na — Nx = M. By using the expression 132(1 . we obtain from (|27|l and l(28|l the equations 



1 dmi (t) 
a dt 



M 2 ~m 2 (t) 



and 



1 dm,2{t) 



a 



dt 



= M 2 - m 2 (t) [1 + 6toi(*)] + 2mi (t)[M 2 + 2m\{t)\, 



As expected, one can easily show that the stationary stable solution is = M. The time dependence has 

been calculated numerically. FIG. [S] shows the dependence of the difference between the exact and the "normal" 
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FIG. 5: Dependence of the difference between the exact and approximate solutions on u = at. The approximate solution corresponds to 
the closure ft3(t) = 0. 



approximate solutions on u — at. One can observe that the difference is surprisingly small even if few particles 
(M = 5) are participating in the process, however one cannot conclude that the closure K^(t) = is always the best. 



IV. LIFETIME OF THE SYSTEM 



Assume that the system S is "live " when it contains at least 1 particle of type X and the number of particles A is 
larger than zero. In the course of the process the X particles are converting the A particles into X particles until A 
particles can be found in the system. Denote by 9 that random moment when the last A particle is converted into 
X. The random variable is called lifetime of the system S. It is obvious that the event {0 < i] is equivalent to the 
event {£(£) = A^ISo}, where Na and Nx are larger than zero. Consequently 

V{9 < t\S }} = W(t\S ) = V{i{t) = N A \S } =p(u,N A \S ) (33) 

is the probability that the lifetime of the system S is not larger than t provided that the initial state of the system was 
Sq = {Na, Nx}- We can see that the properties of the system lifetime are determined by the probability p(t, Na). 3 
The probability density function of the lifetime is nothing else than 

dp(t,N A ) 

w(t\So) = j t , (34) 

and we can see that 

POO 

/ e~ st w(t\S ) dt = w(s\S ) = s U Na (s), (35) 
Jo 

since p(0, Na) = 0, if Na ^ 0. By using the formula (|llfl one obtains 



i c > hohi ■ ■ ■ hx A -i fQa , 

(s + h ){s + An) • • • [s + Hna-i) 



and it is evident that 



oc 



w(0\S ) = / w(t\S ) dt = 1. 



o 



3 Here we omitted the notation referring to the initial state. 
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FIG. 6: Probability density function of the lifetime for three initial states: {10, 1}, {10,5} and {10, 10}. 



It can be worthwhile to study the temporal behavior of systems containing at t = only 1 particle of type X and 
Na > 1 of type A, i.e. to follow the process from the initial state So = {Na, 1} to the dead state Sd = {0, Na + 1}- 
Since in this case Na — n c + 1 and Nx = 1, if n c — 2fc c , where k c is a nonnegative integer, then we obtain 



while if n c = 2fc c 



w(s\S ) 
I . then we can write 

w(s\So) = 



hohi ■ ■ ■ hk c 



(s + h ){s 



hohi ■ ■ ■ hk c 



(s + h kc ) 



{s + h )(s + hi) ■ ■ ■ (s + h kc ) 



In order to calculate the probability density function of the lifetime at different initial states, we should determine 
the inverse Laplace transforms of these expressions. In FIG. El we see three probability density functions of the lifetime 
9 belonging to initial states: So = {10,1}, {10,5} and {10,10}. It is remarkable that at fixed initial number of A 
particles the most probable lifetime of the system decreases significantly with increasing initial number of X particles. 

For the sake of determination of the expectation value E{#} and the variance D 2 {6*} of the system lifetime 9 define 
the function 



K(s\S ) = log w(s\S ) 



N A -1 

E 

71 = 



log 



s + h„ 



We obtain immediately that 



and 



E{9\N A ,N X } = - 



B 2 {9\N A ,N X } = 



dK{s\N A ,N x ) 



ds 



N A -1 



s=0 



d 2 K( S \N Al N x ) 
ds 2 



N A -\ 
n=0 *™ 



(37) 



(38) 



(39) 



By using these formulas we calculated the ratios 

E{0\N A ,N x } 
E{6\1,1} ' 



and 



T> 2 {9\N A ,N X } 
D2{0|M} 
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FIG. 7: Dependence of the average lifetime of the system S on the initial number of A particles in the cases of Nx = 1 and Nx = Na- 
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FIG. 8: Dependence of the variance of the lifetime on the initial number of A particles in the cases of Nx = 1 and Nx = Na- 



and plotted them in FIG. and FIG. GO 

At the first moment it is surprising that the increase of the initial number of A particles results in a decrease of the 
expectation value and the variance of the system lifetime. It seems to be true that could be happened in the folktale: 
the "large" disappears faster than the "small" . In fact, the larger is the initial number of A particles the faster is the 
autocatalytic reaction, and this explains that we see in figures. 
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V. CONCLUSIONS 



The main purpose of this paper was to study the stochastic properties of small systems controlled by autocatalytic 
reaction. In order to give exact solution of the problem we assumed the distribution of reacting particles in the system 
volume to be uniform and introduced the notion of the point model of reaction kinetics. In this model the probability 
of a reaction between two particles per unit time is evidently proportional to the product of their actual numbers. For 
the sake of simplicity we used in this paper the notations X and A for the autocatalytic and the substrate particles, 
respectively. 

In order to make the calculations not very complex we have chosen the simplest autocatalytic reaction: A+X — > 2X , 
and constructed a stochastic models for it. The state of the system at a given time moment is completely determined 
by the actual numbers of X and A particles. The system is called living at time t > when the probability of the 
autocatalytic reaction is larger than zero. 

We calculated exactly the probability p(t, n) of finding n new X particles in the system at time moment t, provided 
that at t = the number of new X particles was zero, and the system contained Nx > particles of type X and 
Na > particles of type A. We have shown that the stochastic model results in an equation for the expectation 
value mi(t) of the new X particles which differs strongly from the kinetic rate equation. Consequently, we can state 
that in this case the stochastic model does not support the kinetic law of the mass action. It is to mention here that 
this statement was already published by Renyi Q many years ago. We calculated the difference between the time 
dependencies of mean values corresponding to the stochastic and the deterministic models. 

Moment- closure approximations of two types have been analyzed, and the results of calculations were compared 
with the exact mean values obtained directly from the probabilities p n (t), n = 0, 1, . . . . We found that the "normal" 
approximation which neglects the third and higher order cumulants can be accepted as relatively " good" approxima- 
tion. 

The probability density function of the system lifetime has been also calculated, and it is found that the most 
probable lifetime depends sensitively on the number of X particles being present in the system at time moment t — 0. 

The main conclusion is that the point model of the stochastic reaction kinetics resulted in a deeper understanding 
of the random behavior of small systems governed by autocatalytic processes. 
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